# File name: 	biv_plots.R
# By: 		tw
# Purpose:		3d kernels
# Last Version:	10.02.09

setwd("D:/Nauka/_R")

library(foreign)
library(fBasics)
library(AdMit)
library(rgl)
library(coda)

d <- read.dta("truth.dta")

load("biv_chao_1998_I_final.RData")
load("biv_chao_1998_I_weak_final.RData")
load("biv_hist.RData")
load("biv_chao_1998_I_weak_robust_final.RData")
 
# Variables saved:
# iv.admit,   par.rmit,   iv.admitis,   iv2.admitis,   iv.admitmh
# wiv.admit,  parw.rmit,  wiv.admitis,  wiv2.admitis,  wiv.admitmh
# t.admit,  t.admitis,  t.admitmh,  tw.admit,  tw.admitis,  tw.admitmh

d$Y_IV,D=d$D_heli,Z=d$Z_heli

print(cov(d$Y_IV,d$Z_heli)/cov(d$D_heli,d$Z_heli))
print(cov(d$Y_IV_weak,d$Z_heli)/cov(d$D_heli_weak,d$Z_heli))
print(cov(dat[,1],dat[,3])/cov(dat[,2],dat[,3]))

#densplot(as.mcmc(iv.admitmh$draws[,1]))
#plot(as.mcmc(iv.admitmh$draws[,1]))
#geweke.plot(as.mcmc(iv.admitmh$draws))
#geweke.diag(as.mcmc(iv.admitmh$draws[,1]))
#cumuplot(as.mcmc(iv.admitmh$draws[,1]))

summary(as.mcmc(iv.admitmh$draws))
for (i in 1:4) print(autocorr(as.mcmc(iv.admitmh$draws[,i])))
print(1-iv.admitmh$accept)
q<-iv.admitmh$draws[,1]
print(length(q[q>=0])/50000)



#densplot(as.mcmc(wiv.admitmh$draws[,1]))
#plot(as.mcmc(wiv.admitmh$draws[,1]))
#geweke.plot(as.mcmc(wiv.admitmh$draws))
#geweke.diag(as.mcmc(wiv.admitmh$draws[,1]))
#cumuplot(as.mcmc(wiv.admitmh$draws[,1]))

summary(as.mcmc(wiv.admitmh$draws))
for (i in 1:4) print(autocorr(as.mcmc(wiv.admitmh$draws[,i])))
print(1-wiv.admitmh$accept)
q<-wiv.admitmh$draws[,1]
print(length(q[q>=0])/50000)


summary(as.mcmc(rwiv.admitmh$draws))
for (i in 1:4) print(autocorr(as.mcmc(rwiv.admitmh$draws[,i])))
print(rwiv.admitmh$accept)
q<-rwiv.admitmh$draws[,1]
print(length(q[q>=0])/50000)




#mean(iv.admitmh$draws[,1])
#mean(iv.admitmh$draws[,4])
#sd(iv.admitmh$draws[,1])
#sd(iv.admitmh$draws[,4])
#density(iv.admitmh$draws[,1])
#density(iv.admitmh$draws[,4])


'HIST' <- function(v1,v2,N){
	l <- length(v1)
	s1 <- mean(v1)-3*sd(v1)
	e1 <- mean(v1)+3*sd(v1)

	s2 <- mean(v2)-3*sd(v2)
	e2 <- mean(v2)+3*sd(v2)

	bp <- matrix(0,N,N)
	
	g <- matrix(0,N+1,2)
	g[1,] <- c(-Inf,-Inf)
	g[2:N,1] <- seq(from=s1,to=e1,by=(e1-s1)/(N-2))
	g[2:N,2] <- seq(from=s2,to=e2,by=(e2-s2)/(N-2))
	g[N+1,] <- c(Inf,Inf)

	'CHECK' <- function(vec,s1,s2,e1,e2){
		if (vec[1]>=s1 & vec[1] < e1 & vec[2] >= s2 & vec[2] < e2) r <- 1
		else r <- 0
		return(r)
	}

	v <- matrix(0,l,2)
	v[,1] <- v1
	v[,2] <- v2

	for (i in 1:(N)) {
		for (j in 1:(N)){
			bp[i,j] <- sum(apply(v,1,CHECK,s1=g[i,1],s2=g[j,2],e1=g[i+1,1],e2=g[j+1,2]))*100/l
		}
	}
	results <- new.env()
	results$grid1 <- seq(from=s1,to=e1,by=(e1-s1)/(N-1))
	results$grid2 <- seq(from=s2,to=e2,by=(e2-s2)/(N-1))
	results$var <- bp
	return(as.list(results))
}

t0 <- proc.time()
hist <- HIST(iv.admitmh$draws[,1],iv.admitmh$draws[,4],20)
t1 <- proc.time()
hist.w <- HIST(wiv.admitmh$draws[,1],wiv.admitmh$draws[,4],20)
t2 <- proc.time()
hist.rw <- HIST(rwiv.admitmh$draws[,1],rwiv.admitmh$draws[,4],20)
t3 <- proc.time()
save(hist,hist.w,hist.rw,t0,t1,t2,t3,file="biv_hist.RData")

print(t3-t2,t2-t1,t1-t0)

persp3d(hist$grid1,hist$grid2,hist$var,col="pink",xlab="beta",ylab="pi",zlab="")
persp3d(hist.w$grid1,hist.w$grid2,hist.w$var,col="pink",xlab="beta",ylab="pi",zlab="")
#persp3d(hist.rw$grid1,hist.rw$grid2,hist.rw$var,col="pink",xlab="beta",ylab="pi",zlab="")






#filled.contour(hist.w$grid1,hist.w$grid2,h)

#rgl.bringtotop()
#rgl.snapshot("betapi.png", fmt="png", top=TRUE )
#snapshot3d("betapi.png")

#plot(wiv.admitmh$draws[,1],wiv.admitmh$draws[,4])

#save(hist,hist.w,file="biv_hist.RData")









